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Abstract 

Lattice measurements of spatial distributions of the light quark bilinear densities 
in static mesons allow one to test directly and in detail the wave functions of quark 
models. These distributions are gauge invariant quantities directly related to the 
spatial distribution of wave functions. We make a detailed comparison of the recent 
lattice QCD results with our own relativistic quark models, formulated previously 
for quite different purposes. 

We find a striking agreement not only between our two quark models, but also 
with the lattice QCD data for the ground state in an important range of distances 
r < 4 GeV~^. Moreover the agreement extends to the L = 1 states [j^ = (1/2)"*"]. An 
explanation of several particular features completely at odds with the nonrelativistic 
approximation is provided. A rather direct, somewhat unexpected and of course 
approximate relation between wave functions of certain quark models and QCD has 
been established. 



^Laboratoire de Physique Theorique est une unite mixte de recherche du CNRS, UMR 8627. 



1 Introduction 



The impressive phenomenological success of the quark model can be most directly appre- 
hended in the spectroscopy of hadron states, where a huge number of predicted states is 
actually observed. A well-known example is the rich spectrum of baryons. However, the 
success of quark models is often unfairly restricted to their success in spectroscopy, while 
they successfully predicted many transition matrix elements between hadrons. In some 
cases this basically reflects SU(6) symmetry, like for the ratios of magnetic moments. In 
other cases, instead, the success relies on a detailed knowledge of the wave functions: a 
striking example is that of the Okubo-Zweig-Iizuka allowed decays of radially excited states 
of charmonium which depend on the precise positions of the nodes of radial wave functions. 

Even though the quark model wave functions cannot be related to quantum field the- 
ory in a rigorous way, their phenomenological success requires a more detailed study and 
comparison with quantities that can be computed in field theory. It is our aim to show 
that spatial distributions of bilinear densities, which in quark models are very directly 
connected with wave functions, are indeed in agreement with results obtained in lattice 
QCD when one of the valence quarks is static. For certain models discussed below, we will 
show that this agreement holds true for a wide range of distances. In this way we directly 
demonstrate the relevance of the wave function picture. 

Our study is not the first one in the field. To our knowledge, the first attempts in this 
direction were made in ref. [1] where the authors predicted the form of equal-time Bethe- 
Salpeter amplitudes for heavy-light mesons in the Coulomb gauge. Their findings were 
in good agreement with wave functions deduced from a certain relativistic quark model. 
Although encouraging, the problem with these results was that the quark models do not 
refer to any specific gauge, in contrast to the above lattice calculations. It is true that 
the Coulomb gauge, as used in ref. [1], is often considered as more physical, yet it is not 
compelling to identify it with the quark model treatment. 

Another lattice determination of Bethe-Salpeter amplitudes was reported in ref. [2], 
where a qualitative agreement with quark model expectations was shown for various meson 
states (see their fig. 5 in ref. [2]). Even though they use a Wilson line to preserve gauge 
invariance of the Bethe-Salpeter amplitudes, the relation with quark model wave functions 
is not straightforward. 

A very useful contribution has been provided by the UKQCD-Helsinki group in ref. [3]. 
The idea was to consider a light quark in the field of a static heavy quark, and consider 
spatial distributions of the charge and the matter densities, both computed on the lattice ^ 
and then compared with corresponding predictions made within a certain Dirac-type quark 
model. There are several advantages in considering spatial distributions : 

(1) This directly eliminates the problem of reference to a gauge, since they are gauge 
invariant quantities: one simply calculates the matrix elements of the light quark 
current at a distance |r| from the static quark. 

(2) The same spatial distributions are easily calculable in quark models, and a rather di- 

^In fact, such distributions were already considered by UKQCD in ref. [3] in the computation of the 
axial coupling gB*B-ix in the static heavy quark limit. See also a discussion made in ref. [S]. 
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rect conceptual connection is maintained with the quark model wave functions. This 
connection is especially direct in the case of Dirac-type quark models, since these 
spatial distributions are then simply quadratic combinations of the upper and lower 
components of the Dirac wave functions. This general type of model has precisely 
been considered by the UKQCD-Helsinki group and compared to their measure- 
ments of charge and matter (vector and scalar) density distributions. Other types 
of models, like the Bakamjian-Thomas approach, also yield these quantities rather 
straightforwardly, as we show in the present paper. 

(3) Moreover a connection with physical observables is more direct. A simple integration 
over these densities yields quantities of direct physical interest, such as the pionic 
decay couplings of 1~ 0" (e.g. B* -> Bn), and 1+ 0+ (e.g. B[ -)■ Bqh) often 
referred to as ^ and g [6] . 

We find a few difficulties in the work of ref. [3]. In their lattice calculation, they claim 
to determine an extensive spectrum of static heavy-light mesons, including the higher 
orbital and radial excitations, up to L = 3 and n = 3, with surprisingly massive radially 
excited states. When trying to reproduce these results by quark models, they are led to 
consider rather uncommon confining potentials that are a mixture of scalar and vector 
parts with a very large string tension. With such potentials, they do not reproduce the 
charge distribution within the doublet of lightest mesons, = (l/2)~, which should in 
principle be the one for which the models work best. 

The aim of the present work is to revisit that problem by using more standard quark 
models, based on a scalar linear potential in addition to the Coulomb-like short-distance 
part, a form formulated and used long ago, independently of the present problem. We 
also use the improved lattice computations discussed in the companion paper [7] , including 
a larger spatial extension, better accuracy, and extending the computation to the axial 
distributions for both the ground states [j^ = (1/2)"] and the nearest orbital excitations 
[J^ = (1/2)1- ' 

As we will show, the two types of models chosen in this paper give similar results, and 
are well compatible with the lattice data. In particular the (vector) charge distribution for 
(l/2)~ states is now very well reproduced, while the comparison between the lattice results 
and the distributions deduced from the quark models is more than satisfactory not only for 
the (l/2)~ doublet, but also for their excited counterparts with = (1/2)"''. We should 
also emphasize that in this paper we restrict our attention to elastic matrix elements: we 
do not consider distributions of the transition matrix elements between states belonging 
to different doublets of heavy-light mesons, or between radially excited and ground states 
(for the spectrum of such states please see ref. [S] ^). 

^Throughout this paper we will use the heavy quark spectroscopic notation = (1/2)" which stands 
for the ground state doublet of heavy- light mesons = [0~,1~], while the lowest excited states with 
L = 1, correspond to = (1/2)+, i.e. to the doublet of heavy-light mesons = [0+, 1+]. 

■^Thc authors of that paper also plan to compute the transition matrix elements of the axial charge 
between radially excited and ground heavy-light static mesons. 
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2 Two proposed quark models 



We will confront the lattice data with results obtained by using a Dirac-type quark model, 
and one based on the Bakamjian-Thomas (BT) approach [9]. We advocate them both, 
because both exhibit a large set of interesting features and successful results. Although 
their respective structure seems rather different, they share several common features. 

First, both quark models are relativistic, i.e. they have a relativistic kinetic energy. By 
now, there is no doubt that in contrast to heavy quarkonia, heavy-light systems are highly 
relativistic, necessitating at least a partially relativistic treatment of kinetic energy p!0] . 
In the two quark models we consider here, this is made by (a) using the standard Dirac 
equation, or by (b) introducing the free kinetic energy of a quark, y^p^ + m^, as in the 
well-known potential model by Godfrey and Isgur (GI) [11]. Other relativistic features 
come from the treatment of the potential as having a Lorentz structure of vector or scalar 
exchange in the space of Dirac spinors. 

Second, they are based on the standard assumption of a linear confining (scalar) poten- 
tial, plus a Coulomb-like short- distance part. We are of course aware that many objections 
have been formulated against a scalar confinement potential. Phenomenologically, how- 
ever, it offers a simple and satisfactory description of the data, for which there exists no 
real alternative. ^ 

Other than these common features, the two models are rather different, and so are the 
expressions for hadronic matrix elements, which makes the agreement of the physical results 
obtained from both quite surprising (see sec. 13.31 below). Before entering the specifics of 
both models, we should add a few comments. 

(i) The meaning of the Coulomb-like potential: As usual, a Coulomb-like potential is in- 
troduced to mimic physics at "short-distances". It is, however, unclear what that 
assumption exactly means. It should be stressed that the naive idea of one-gluon 
exchange dominating at short-distance is not a sufficient justification for the assump- 
tion. Indeed, studies of the static potential by using lattice QCD [H] suggest that the 
one-gluon exchange can dominate only at very small distances, hardly accessible from 
the lattice data, and only a small part of the static potential is describable by pertur- 
bative QCD, even after including the perturbative higher order corrections. On the 
whole, the Coulomb-like potential in phenomenological models covers a large range 
of distances and should be probably considered only as a phenomenological descrip- 
tion that, together with the linear contribution, provides the medium-range potential 
responsible for the bound states. This is manifested by the fact that k, = 4as/3 is not 
necessarily small. In the Dirac model, as well as in the Cornell model of quarkonia, 
it is in fact large. 

(ii) Interest of the Dirac-type models: We may say that the Dirac equation formalism 

^Recently, a new phenomenological objection has been formulated by the Wisconsin group [12]. They 
predict a rapid variation of the excitation energy with the light quark mass in a heavy-light system, which 
would contradict the experimental data. We find this to be a consequence of the fact that they assume 
a hard Coulomb potential, cx 1/r. In fact their prediction changes completely after softening this 1/r 
behavior, as should be done according to asymptotic freedom. As a result one finds that the excitation 
energy is slowly varying with the light quark mass, as observed. 
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has two main features which fit particularly well the present problem. First, it 
corresponds to the static limit, which is the one considered on the lattice side. Second, 
the expressions for various spatial distributions of densities is particularly simple; they 
are just quadratic combinations of the upper and lower components of the radial Dirac 
wave functions. A phenomenological drawback of the static limit is related to the 
fact that the spectrum of states is not the physical one, and one cannot directly 
rely on the physical spectrum to fix the model parameters. One instead has to rely 
on extrapolations from the sectors with a heavy c- or 6-quark to the static limit. 
Moreover, the experimental checks of the predictions made in the static heavy quark 
limit can be at best made in these heavy-light systems. 

(iii) Interest of the Bakamjian-Thomas approach: This approach allows one to formulate 
relativistic states in motion, for any mass of all the quarks. In the BT-type models 
one has an exact representation of the Poincare group. However, it is not known 
in this scheme how to define fully covariant current operators. Interestingly, it was 
found that, in the heavy quark mass limit, and the current acting on the heavy quark, 
covariance of current matrix elements is obtained, as also Isgur-Wise scaling [28]. 
In the case of finite mass, and in the present case of density distributions, where the 
active quark is the light quark, full covariance is lost. However, we hope nevertheless 
that important relativistic effects are still incorporated in the model. The BT scheme 
considerably enlarges the domain of physical checks that can be made, although here 
we work in the static limit. In addition, we adopt as the rest frame wave equation 
the relativistic wave equation derived by GI [H], which in our opinion is unique in 
having been tested over a particularly large set of states and reactions. 

Concerning the parameters of both the Dirac and BT models, we must first say 
that neither of them has been fitted from the spatial distributions themselves, that 
constitute predictions of the two schemes. In both models the parameters have been 
fixed to describe the hadron mass spectrum only. In the case of the BT model we 
adopt the original parameters of the Godfrey and Isgur (GI) potential 

2.1 A model with Dirac equation 

Our Dirac model is, of course, based on the equation [17] 



The differences between the eigenvalues E correspond to the differences between the 
hadron masses, and are used to fix the parameters of the potential V{r), which is dis- 
cussed below. Once these parameters are fixed, we only deal with eigenfunctions \I^(r). 
From the knowledge of the quadrispinor \l/(r), one computes straightforwardly the density 
distributions, which are defined by 



a ■ p + m/3 + V{f) 



E ^(f) = 0. 



(1) 



P0(f) = ^(f)1"C^(f) , 
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with the desired Dirac matrix being O = l,a,/3, for vector [V), axial (A), and scalar (S) 
density, respectively. One must also take into account the angular momentum j of the light 
quark and of the heavy quark spin, to form the desired - eigenstates. It is then very 
easy to establish, by purely algebraic calculation, the explicit expressions for the various 
matrix elements of local densities for the ground states, in terms of radial wave functions 
f[/\ir) and g^L{r), which enter the full wave function as 



fHr) yUr) 



(3) 



Notice that yj'mir) denotes the combination of spherical harmonics with spinors standardly 
used in the solution of the Dirac equation, with the definitions of Landau and Lifshitz |16] : 
k is the Dirac quantum number [A; = / for / = j + 1 , and k = — (/ + 1) for / = J — |] , which is 
k = —1 {k = 1) for the = 0~ (0^) state, and k = —2 for = 2"*". After inserting this 
into the Dirac equation, and integrating over the angular variables, one arrives at a system 
of coupled differential equations for the radial functions fj{r), gj{r), that can be computed 
numerically by a method of expansion over a basis of orthogonal functions explicitly given 
in the Appendix of the present paper. In what follows, the results are obtained in the 
basis of dimension n = 40. ^ The details of such a calculation have been already described 
in detail in ref. [17] where we obtained the expression for the matrix element of the axial 
current. Here we give the results for all three matrix elements: 

Pv{r) = l/r/aWT + kr/a^^^r ' 
PA{r) = \A~/l{r)\'-l\9;/\{r)\\ (4) 
Psir) = \fylir)\'^ - \gy2{r)\'^, 

for the vector, axial, and scalar density distributions, respectively. ^ As mentioned above, 
these expressions are simple quadratic combinations of the radial wave functions. One must 
warn that, in the case of the axial current operator a, since it is a space-vector, one must 
consider the transition between 1~ and 0~ states. One could define various distributions 
according to the direction of a. We choose the one calculated in the companion paper [7], 
which corresponds to the emission or absorption of the pion between the 0~ and the 1~. 
Moreover, there are sign conventions to be fixed (emission or absorption, relative phase of 
the two states). The above one is in agreement with what has been chosen in ref. [7]. It is 
trivial to see that there is a linear relation among the above three distributions, namely, 

Ps{r) = ^ [3pA(r) - py(r)] , (5) 

"'In the Appendix we provide the reader willing to reproduce our results with expansion coefficients 
in the basis of dimension n ~ 15. Extension to n = 40 is obviously straightforward. If necessary the 
coefficient in the basis with n — AO can be obtained upon request from the authors. 

^Note that the vector and scalar densities are often referred to as charge and matter densities, respec- 
tively. 
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that is also valid in the Bakamjian-Thomas approach, as we show in the next subsection. 
This feature is manifest in models but not in QCD where the scalar density is renormal- 
ization scheme and scale dependent, in contrast to the vector and axial ones. With that 
in mind we will also consider the scalar density distribution and check the extent to which 
the above relation is valid (up to a proportionality constant). 

Similar quantities for the = (1/2)+ doublet have been computed on the lattice 
too [7]. The corresponding quark model expressions are: 




where the superscript + denotes parity and is used to distinguish from the expressions 
given in eq. (jl]). ^ Like in the case of ground states, also these densities in quark modes 
satisfy a simple relation, 

pUr) = l[pUr)-3p\{r)] • (7) 

To obtain the above radial expression in the axial case requires consideration of the tran- 
sition between l"*" and 0"'' heavy-light mesons, once more because the operator is a space- 
vector. The convention chosen here for the axial distribution is opposite to the one in the 
lattice paper. We take this into account in our comparison below, by converting the sign 
of the lattice data. 

At this point we would like to underline the fact that the vector density is of course 
positive definite, while the axial and scalar densities are not. In fact they are expressed in 
terms of differences between the large and small component contributions [see eqs.dH |6])], 
and therefore can change sign as a function of r. This actually happens for = (1/2)''" 
states as we show below, and clearly demonstrates the crucial importance of the relativistic 
treatment since in the nonrelativistic approximation the small component vanishes. See 
Fig. I2.1.1l where the three densities obtained from the Dirac model are compared with those 
obtained from lattice QCD. 

One should be particularly careful in choosing the orthogonal basis to diagonalize the 
Hamiltonian when dealing with (1/2)+ states. It was proposed in ref. [IB] to use the 
pseudo-Coulombic (PC) basis, and useful formulas were provided. These formulas were 
also used in ref. [17]. Near the origin the basis functions behave as ~ r^pc ^ However, in 
[T8] it was assumed that the orbital angular momentum ipc of the basis is identical to 
the i of the Dirac equation {i = j ± 1/2) which gives the correct behavior for the upper 
components, with a regular potential. However, the lower components behave differently, 
as , where i' = i^l. In the case L = 1, j = 1/2, the i of the Dirac equation i = L = 1, 
but i' = 0, and therefore the ipc of the basis functions for the lower components should 
be now chosen to be = 0. This has been implemented in the present paper. The main 

^To avoid ambiguities, when we discuss the features apphcablc to both the ground and excited states, 
we wiU denote them by a superscript (+), e.g. Py'\r) 
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effect is that the above densities do not tend to zero, having always a contribution from 
the lower components. Our previous results reported in ref. [17] for the spectrum and the 
matrix elements remain essentially unchanged. 

2.1.1 Parameters of the Dirac model 

The remaining problem is then to determine the potential and the mass in eq. ([H). We 
choose the following form of the potential: 

V{r) = p{ar + c)- , with d{r) = 1 - g-^"/^)' . (8) 

The function d{r) in the vector part of the potential is chosen to temper the Coulombic 
singularity at the origin, and its introduction does not alter the discussion made in our 
previous paper in which we took d{r) = 1 [17]. As for the scalar part of the potential its 
form is the same as in ref. [17], where we discussed the puzzling discrepancy between the 
QCD sum rule prediction of ^ and its much larger experimental value, 

9 = 9a = Jd^fpAir). (9) 

Notice in particular that the constant c in eq. ([H]) combines with the constituent quark 
mass into a constant term which does not have to be positive. Since our paper [T7|, new 
information concerning the spectrum of heavy-light mesons has been obtained both in 
experiment and through numerical simulations of QCD on the lattice. This is why we will 
somewhat readjust our set of preferred parameters given in ref. [17]. The most important 
constraint on the model parameters comes from AEi = mi^L=i) — mi^L=o) ~ 0.4 GeV. 
Because we are working in the heavy quark limit, we should in principle distinguish the 
(1/2)"'' states from the (3/2)''' ones, both being L = 1. The mass degeneracy of these 
states in the heavy quark limit is the second important constraint on the choice of model 
parameters. We therefore require the spin-orbit splitting to be negligible, in agreement 
with the celebrated model by Godfrey and Isgur [11] . In the Dirac model this implies, with 
the known value of the string tension, a, that the Coulombic term is rather large, far from 
any perturbative regime. This is not a problem, in our opinion, since we do not consider 
it as a short-distance one-gluon term, but rather as an effective contribution at medium 
distances. 

As already mentioned, only the combination 

Co = m + c, (10) 

enters the Dirac equation and it may be negative. From experience with potential models, 
we know that the confining potential is most often written as ar + c, with c being negative 
and large. Therefore, one should not object to cq < 0, as it results from our fits. In ref. [T7| 
we found cq = —0.2 GeV, which was an important distinctive feature of our model. As we 
explain below, in this paper we find cq = —0.1 GeV, still negative. The fact that cq < 
also implies an important cancellation of the linear term of the potential in the region 
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r [GeV"'] r [GeV"'] 

Figure 1: Radial distributions for the light-light vector and axial current, as well as for the scalar density, 
in the case in which the external states belong to the lowest lying doublet of static heavy-light mesons 
[j^ = (1/2)^] is presented in the left plot, while in the right one the external states are the static heavy- 
light mesons belonging to the nearest orbital excitations [j^ = (1/2)"'']. The results are obtained by using 
the Dirac-type model discussed in the text. 



around r ^ 1 GeV~^, for a ^ 0.2 GeV^, which will play an important role in the apparent 
approximate 75-symmetry discussed in sec. 13.31 

In the following we discuss in more detail the data used to determine the parameters 
in the potential ([8]). 

• Experiment: The experimental observations of D and B states with L = 1 [j = 
(3/2)"''] [ini [20] allow us to linearly extrapolate to the static heavy quark mass limit 
where we find AEi ~ 0.4 GeV, for both strange and nonstrange heavy-light mesons. 

• Lattice QCD: A recent lattice QCD study of the spectrum of heavy-light mesons in 
the static limit reports AEi ~ 0.4 GeV, and indicates the spin-orbit splitting to be 
negligible, both in quenched approximation and after including A^^ = 2 dynamical 
light quarks [21]. ^ 

Notice, however, that our assumption of small spin-orbit splitting has been challenged in 
the Dg sector: m^,*^ — ttlu*^ ~ 0.25 GeV. ^ Presumably, what is observed experimentally 

^Note that the conclusions of the ETM Collaboration are different. Their first paper [22] is superseded 
by the second one [23], which interprets differently the strange state, by a "partially" quenched calculation 
where the sea is constituted hy u/d loops. In ref. [23], the spin-orbit mass difference is the same for strange 
and nonstrange states, i.e. 0.1 GeV, thus larger than the negligible one found in ref. [21]. Note that the 
calculations of ref. [H] are unitary" , i.e. the sea and the valence light quarks are degenerate. In the case 
of the non strange mesons, the ETMC data are also unitary, so the two results are comparable but do not 
agree. We have no explanation for this discrepancy. 

^Note that such a large value is also found in the case of non-strange D mesons: m^i* — m^j. ~ 
0.14(3) GeV. 
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for the Dgj states, i.e. very large spin-orbit splittings, is a finite mass effect, although that 
phenomenon is not easy to understand in quark models. 

Another important spectroscopic information is a newly identified radial excitation of 
Dg [21]. The radial excitation energy has been determined cleanly on the lattice too. In 
the static limit, in quenched approximation and for the light quark around the strange 
quark mass, it is found to be around A£" ~ 0.6 GeV [8], in qualitative agreement with 
experiment and with quark models for the Dg. Similar methods applied to Np = and 
Np = 2 lattice data [21] give a similar value. In summary, the radial excitation is found 
to be higher than the L = 1 states. 

We emphasize that the parameters used in this paper are not far from the set of pa- 
rameters used in ref. [17]. We have in fact tried to remain as close as possible to this old 
set of parameters, while maintaining the assumption of a very small spin-orbit splitting, 
which leads us to the following modifications: 

(i) To increase AEi while maintaining the spin-orbit splitting small requires a simulta- 
neous increase of the string tension a, and of the strength of the Coulombic term k, with 
respect to the values we obtained in ref. [17]. We also choose to respect the mass differ- 
ence between the strange and nonstrange heavy-light mesons, ~ 0.100 GeV, to obtain: 
a = 0.25 GeV^, k = 0.75, with d{r) ^ 1 only at very small distance. The constant cq is 
found to be about cq = —0.1 GeV for the strange, and cq = —0.35 GeV for the nonstrange 
states. 

(ii) The choice of function d{r), introduced in eq. ([H]), is meant to maintain a small 
splitting between (3/2)"'" and (1/2)"*" states, independently of the light quark mass. At small 
distances d{r) <^ 1 allows one to kill the Coulomb singularity. With a strict Coulomb-like 
potential, and large n, we would have obtained the spin-orbit splitting increasing too rapidly 
with the light quark mass, which contradicts the observations made in lattice QCD [2Tj . 
Indeed the larger the light quark mass, the more the wave function is concentrated at the 
center and the weight of the Coulomb piece in the potential increases. As a consequence 
the contribution of the vectorlike part to spin-orbit splitting exceeds the scalarlike one, and 
the spin-orbit splitting grows with the light quark mass. To suppress that effect, it suffices 
to cut off the singularity, which we do by introducing a Gaussian cutoff, namely a factor 
d{r) = 1 — exp[— (r/i?)^]. We choose R ~ 0.25 GeV~^, a rather short-distance. This cut is 
enough to keep the spin-orbit splitting very small independently of the light quark mass. 
A similar procedure to soften the Coulomb singularity has been used in the GI model. 

With the above set of parameters and for the strange light quark we obtain AEi ~ 
0.4 GeV and for the spin-orbit splitting Ep^^^ — Ep^^^ = —0.01 GeV, as expected. We then 
deduce the mass difference between the radial excitation and the ground state to be AE' = 
0.576 GeV, in agreement with lattice and experimental findings. For nonstrange states we 
have AEi ~ 0.38 GeV and Ep^^^ - Ep^^^ = 0.01 GeV, and predict AE' = 0.556 GeV. The 
difference between the strange and nonstrange ground states is 0.1 GeV. 

^°It is important to mention that the conversion from lattice to physical units in both ref. ^ and ref. [21] 
has been made by using rg = 0.5 fm. 
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2.2 The Bakamjian— Thomas approach to relativistic quark mod- 
els 

In this Section we compute the same distributions of heavy-hght mesons but in the class 
of Bakamjian-Thomas (BT) relativistic quark models [U [251 ISSl |271 128]. This is a class 
of models with a fixed number of constituents where the states are covariant under the 
Poincare group. The model relies on an appropriate Lorentz boost of eigenfunctions of a 
Hamiltonian describing the spectrum at rest. As it has been explicitly shown, the model 
has the important properties that the matrix elements of currents between hadrons are 
covariant in the heavy quark limit and exhibit Isgur-Wise scaling [281 ISHl |30l |31] . Given 
a Hamiltonian describing the spectrum, the model provides an unambiguous result for 
the elastic Isgur-Wise function C,{w) [2H1 132], and satisfies several sum rules (SR) in the 
heavy quark limit of QCD, such as Bjorken or Uraltsev SR, as well as SR involving higher 
derivatives of ^{w) at zero recoil [331 [M]. 

One can compute the V, S, and A distributions for mesons within the BT scheme using 
a Hamiltonian describing the spectrum, like the GI model [H], that has been used else- 
where [29] to compute the elastic Isgur-Wise function ^{w) and the inelastic ones ti/2{w), 
T'3/2{w) [32] • Since the computation of the light-light matrix elements within the BT scheme 
has not been presented before, we present it here in more detail. 



2.2.1 The Bakamjian-Thomas approach to quark models 

We first recall some basic elements of the BT-approach to quark models. The con- 
struction of the BT wave function in motion involves a unitary transformation that relates 
the wave function \l/if^L,s„(pi, ■ ■ ■ ,Pn) in terms of one-particle variables, the spin Si and 
momenta Pi to the so-called internal wave function ^^(-P, k2, - ■ ■ , k^) given in terms 

of another set of variables, the total momentum P and the internal momenta ki, k2, ■ ■ ■ , kn 
(X] — 0) PH]- This property ensures that, starting from an orthonormal set of inter- 

i 

nal wave functions, one gets an orthonormal set of wave functions in any frame. The 
basis ■ ■ ■ ,pn) is useful to compute one-particle matrix elements like current 

one-quark matrix elements, while the second s„(-P) ^2, ■ ■ ■ ? kn) allows one to exhibit 

Poincare covariance. In order to satisfy the Poincare commutators, the unique requirement 
is that the mass operator M, i.e. the Hamiltonian describing the spectrum at rest, should 
depend only on the internal variables and be rotationally invariant, i.e. M must commute 
with P, ^ and S. The internal wave function at rest (27r)^(5(P)y9s^^...^s„(^25 ■ " " ? ^n) is 
an eigenstate of M, P (with P = 0), S*^ and Sz, while the wave function in motion of 

momentum P is obtained by applying the boost Bp, with P = + involving the 
dynamical operator M. 

The final output of the formalism that gives the total wave function in motion, 
in terms of the internal wave function at rest (psi,--- ,s„ (^2? ■ ■ ■ ? kn) is encoded in the following 
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formula 




X Yl [A(Ri)],^,,/<^4,..,s;(A^2,--- ,4), (11) 

where p° = y^pf + mf, and Mq is the free mass operator, given by Mq = ^(^PiY- 

The internal momenta of the hadron at rest are given in terms of the momenta of 
the hadron in motion by the free boost ki = B^^ pi where the operator Bp is the boost 

i 

(a/p^, 0) — )■ p, the Wigner rotations Rj in the preceding expression Rj = Bp.^B^^p_Bfc. and 

i 

the states are normalized according to (P', S'z\P, Sz) = (27r)^5(P' — P)6s^,si,- 

The current one-quark matrix element acting on quark 1 between two hadrons is then 
given by the expression 

- ,Pn)* <p'i,s'i|/')|pi,si > M/f^,...,,jpi,--- (12) 

where \l/f^ ... ^^(pi, ■ ■ ■ ,Pn) is expressed via the internal wave function as in eq. ( fTTl) . and 
(p'l, s'^l J^^^Ipi, Si) is the one-quark current matrix element. 

As demonstrated in refs. [28l [30] . in this formalism, current matrix elements in the 
heavy quark limit are covariant, they exhibit Isgur-Wise scaling, and one can compute 
Isgur-Wise functions like ^{w), ti/2{w), T3/2{w) [29] . 

In the present paper, we deal with light quark currents between heavy-light meson 
states, and there is no a priori reason for the result to be covariant. To compute the 
desired radial distributions we need to specify the mass operator M, which we choose to 
be the one appearing in the GI model, and we work in the meson rest frame. 



2.2.2 Ground state vector, scalar, and axial densities 

We need here to compute the matrix elements of the operators O = qVq where F is 
any Dirac matrix, in practice F = 7*^, 1, 77^, for the vector current, scalar, and axial vector 
current densities respectively. In ref. [28] we did compute matrix operators of the type 
f ll2p where we took the convention "1" for the heavy and "2" for the light quark. The 
heavy mass limit was taken for the quark 1. In the case we are treating here, the light 
quark is coupled to the current, while the heavy mass limit is to be taken for the spectator. 
We keep the same notation as in ref. [28], and for finite mass for both quarks the matrix 
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element f lT2|) one writes 



(i?(P')|g(0)rg(0)|i?(P)) 
dpi 1 



(27r)3p0 



1,01,0 



^2^2° V (^1 + '^Ol^" + ^2) V (^'1 + mi){k''^ + 1712) 



-Tr [r (m2 +^2) (1 +^) ("^1 +A) (1 +/) K > 
where the internal wave function at rest, v^(/c2), is normalized as 

dk2 



(13) 



(2vr)= 



(14) 



Pi, p']^ and p2i p'2 ("^1 s-iid 7712) are the quark four- momenta (masses) of the heavy and light 
quarks respectively. The four-vector u, Mq, and the momenta ki and pi are related as 



u 



P1+P2 
Mo 



Mo = ViPi+P2y 



Buki =Pi {i = 1, 2) 



(15) 



where ki and /c2 are the four-momenta of the quarks in the rest frame of the B meson 
and Bu is the boost associated with the four-vector u. Taking the heavy quark limit 
(mi — )■ 00) in eq. (fT3!) amounts to the following replacements: Mq — )■ mi — )■ m^, m*-'^ — ?■ w*-'^ 
, pf^ /mi — t^^'^ ki^^ /mi — > 1, /cg'*" — ^ P2'* ' '^'''^ where w and f' are the initial and final 
S-meson four- velocities. In the frame in which P = 0, P' = q one then obtains: 
(i) Vector current, 



{BmmiMomo)) 



dk 
(2^ 



V{k + q)* ip{k) 



(^/W+^m? + m){\J (A; + g)2 + m^ + m) + k- {k + q) 



\/A;2 + m^-y/ (k + q)'^ + m'^{ + + m) ( \J {k + q)"^ + m^ + m) 



(16) 



where we write m2 = m, since the light quark mass is the only mass in the game. The above 
matrix elements satisfies (i?(0)|g(0)7''g(0)|i?(0)) = 1, as expected from gauge invariance. 
(ii) Scalar density, 

dk 



(^/W+^m? + m){\/ [k + g)^ + + m) — k ■ {k + q) 



+ m'^sj (k + q)"^ + m^('\/p + m^ + ^){\J {k + g)^ + m^ + m) 



1 ) 
2 



(17) 



^-'^As noted above, the BT states are normalized acccording to {Bd{P')\Bd{P))'BT = (27r)'^(5(P — P'). 
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which leads to the scalar coupling, 

9s = (BiOMOUOmO)) ■ (18) 

The Fourier transforms of f|T6l) and f|T7|l finally give the desired densities in configuration 
space 



Pv[r) 



dq 

(M 



dq 



.e^«"(5(g-)|g(0)g(0)|5(0)).(19) 



The densities depend only on the modulus r = |r| and, after a lengthy integration over the 
angles, one gets 

Pv{r) = pi(r) + p2{r) , ps{r) = pi(r) - p2{r) , 



where 



2^2712 r 7o 



dq q (p{q) 



sin(gr)( y/gM-rn^ + m) 



\ 



2^/271^ Jo 



dq q if{q) 



\J q^ + m?{ q^ ^vn?' + m) 

[sin(gr) — qr cos(gr)] 

m] 



(20) 



(21) 



(22) 



y ' + ( -\/g2 + + m) 

both square-functions of r, which is expected because f lT6|l and f lT7|) are written as a linear 
combination of two terms, convolution products of the same functions in momentum space, 
(iii) Proceeding in the same way, for the matrix element of the axial current 

(fi(g-)|g(0)(7-n)7'g(0)|r(0,6-)), (23) 

one finds two rotationally invariant terms pA^f) = Pa''('")(^ ■ ^ + Pa''('")(^ ' ^){^' fo^' 
which we obtain simple results, related to P5(r) and pvi'f') from eq. (1201) . 

(^) = Psir) , piJ^(r) = pvir) - ps^r) . 
In configuration space one recovers the normalizations. 



dr Anr'^pyij-) = 1 , 



dr 4nr'^ps{r) = gs . 



(24) 
(25) 



Taking for the axial densities the components such that n = e, one obtains 



PA{r) =pi^^(r) + |p^^^(r) 



.(2)/ 



and therefore one verifies the same relation as in the Dirac model, namely 
and, by integration over r. 



1 2 

PAir) = - pv{r) + - psir) 



9 = -^ + -^9s. 



(26) 
(27) 
(28) 



where we use the notation = Notice that relation (127|) shows that also in the BT 
scheme there are only two independent functions like in the Dirac model dl]). 
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2.2.3 The choice of the mass operator M: The Godfrey-Isgur model 



Let us now particularize the above expressions for the choice of the mass operator M given 
by the GI model [H] that describes the whole set of meson spectra qq and Qg, where g is a 
light quark (g = u, d, s) and Q is a heavy quark {Q = c, b), with the important exception 
of the recently discovered narrow states Dgj O"*" and 1^, that are too low compared with 
the predictions of the model. The model contains a relativistic kinetic term of the form 

K = sjkl + ml + sjkl + ml , (29) 

that is identical to the operator Mq at rest [28], and a complicated interaction term that 
includes : (1) a Coulomb part with a q^ dependent a^, (2) a linear confining piece, and (3) 
terms describing the spin-orbit and spin-spin interactions. All singularities are regularized 
- e.g. terms of the type 5{r) or l/m2, where m2 is the light quark mass. The Hamiltonian 
H depends on a number of parameters that are fitted to describe all the meson spectra. 
By inserting the GI wave functions (see the Appendix of the present paper) in the above 
expressions we computed the densities pv,A,s{''^)i and after integrating over we obtain 

(75 = 0.56, ^=0.71. (30) 

Notice that at large distances r, pAij") becomes larger than ps{r), and since this portion 
of the integral is enhanced by the integration measure growing like r^, the resulting values 
are different. For the mean square radii, defined by 

POO 

(r^), =An dr r^p,{r) (z = V, S, A) , (31) 

we obtain 

{r^)y = 5.174 GeV^^ {r'^)s = 1.721 GeV"^ (r^)^ = 2.872 GeV~^, (32) 
that satisfy relation ( 1271) . 



2.2.4 Vector, scalar, and axial current densities for the lowest excited states 

(1/2)+ (L = 1) 

Like in the discussion with the Dirac-type quark model, we now extend our discus- 
sion to the excited heavy-light static mesons with L = 1, and consider the corresponding 
transitions (1/2)+ (1/2)+, i.e. 0+ 0+ and 1+ 0+. 

After proceeding like in the previous subsection, we obtain the following expressions : 



(i?(0+)(g-)|g(0)7°g(0)|5(0+)(0)) 



dk V*i/2ik + q) ipx/2{k) 



(2vr) 



\k\ \k + q\ 



n 1/2 



^2 _|_ _|_ rn){y {k + q)"^ + + m) 



+ m? \ {k + q)"^ + m? 



( Y /c^ + m? — m){y {k + q)'^ + m? — m) + k.{k + q) 



(33) 
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(S(0+)(g-)|g(0)g(0)|B(0+)(0)) 



dk ^*i/2ik + q) ^1/2$) 



(2vr) 



\k\ \k + q\ 



k'^ + m? + m){\j (/c + g)2 _|_ ^2 _|_ 



1/2 



X 



(V fc^ + m? — m){y {k + q)'^ + m? — m) — k.{k + q) 



which, when Fourier transformed, lead to 

Pvi^) = Ptir) + ptir) , pj(r) = -p|(r) + p+(r) 
with Pi{r) and p^(r) given by the squares. 



(34) 



(35) 



p+(r) 



P2 I'^J 



1 1 
2V^r Jo 



dq <-Pi/2{q) sin(gr)(v/g2 _^ ^2 



V^i/2(Q')[sin(gr) — grcos(gr)] 



\/ q^ -\- vn? + m 



(36) 



As far as the axial current is concerned, by Fourier transforming the matrix element de- 
scribing the 1+ 0+ transition, (S(0+)(g)|g(0)(7 ■ n)7^g(0)|5*(l+)(0, e)), we find 



Pa(^) = Pa ^^\r){n ■ e) + ^'^\r){n.f){e ■ r) , 
but with relations to ^^(r) and Py(r), that are different from eq. i^^, namely. 

Pa ^'\r) = -pUr) , p1 ^'V) = P^r) + P^(r) , (3. 
that finally yield the result for the couplings 



(37) 



9 3 ^9s 



(39) 



in the notation frequently encountered in the literature, 'g = g\. As in the previous sub- 
section, for the numerical calculations we use the GI potential model and, after integrating 
over r, we obtain 

r7+~0.30, 5^ =(7+ -0.13. (40) 

On the other hand, for the mean square radii, defined in (13T1) . we get 

(r2)+ = 5.42 GeV-2, {r^s = 2-73 GeV-^, {t^a = - 0-015 GeV-^ , (41) 

that, of course, satisfy a relation of the type ( 139|) . All three densities obtained by using 
the BT model are compared with the lattice QCD results in Fig. 12.2.41 
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r [GeV"'] r [GeV"'] 



Figure 2: Plots similar to those presented in fig. \2.1.1\ but obtained by using the Bakamjian-Thomas 
type quark model, discussed in the text. 



3 Comparison of the predictions of the quark models 
with the lattice spatial distributions 

3.1 Short description of the lattice results 

In this section we briefly remind the reader of several details concerning the recent results 
for the radial distributions discussed in this paper as obtained from the QCD simulations 
on the lattice. For more details, please see ref. [7]. The calculation of spatial distributions 
is performed with a static heavy quark source, and with the light Wilson improved quark 
action. The gauge field configurations, containing Np = 2 mass degenerate light quark 
flavors, were obtained by using the Iwasaki-type of the lattice Yang-Mills action [15]. The 
bare coupling constant go is chosen according to (3 = G/qq = 2.1, and the associated value 
of the lattice spacing is found to be around ^ 2 GeV. Although dynamical quarks 
are considered, no attempt has been made to include disconnected diagrams. To improve 
the signal-to-noise ratio a suitable hyperblocking of the heavy quark propagator (Wilson 
line) has been implemented. The lattice size used in that computation, 24^ x 48, means 
that the spatial extension of the lattice is 12 GeV~^, so that due to periodic boundary 
conditions, the points up to radial distances 12 a ~ 6 GeV~^ could actually be explored. 
Also, all directions for the spatial points of the lattice have been exploited, allowing for a 
large averaging and giving a very dense set of points except at small distance. The four 
values of the bare light quark mass (each distinguished by its Wilson parameter Kg) have 
been chosen around the strange quark mass. Moreover the valence and the sea-quark 
masses have been kept degenerate, i.e. the quarks in loops have consistently the same mass 
as the valence ones. The physical densities are obtained from the ones computed on the 

i^More specifically Kg e {0.1357,0.1367,0.1374,0.1382}, with Kg = 0.1374 being the closest to the 
strange quark mass. The critical Kc = 0.138984(13) corresponds to the chiral limit [T5] . 
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lattice via 



Po{r) = l/a^ ZoPTir/a), 



(42) 



where a is the lattice unit (lattice spacing), pQ**(r/a) stand for the bare lattice data with 
O e 5", A, and Zoig^) is a relevant renormalization constant, all being determined in the 
same lattice calculations. The renormalization constants Zv^AiOo) have been calculated 
nonperturbatively through Ward identities, and for = 0.1374, their (9(a)-improved 
values are Zy = 0.814, Za = 0.847. Since these values were found to be consistent with 
what one obtains by using the boosted perturbative estimates, we use the perturbatively 
evaluated Z^^{2 GeV) = 0.639. Notice that the above renormalization constant Zy ensures 
the vector current conservation in the lattice regularization scheme, i.e. 



The corresponding integral concerning the axial current differs from 1, it is often denoted 
as 'g or qa, and has a physical meaning. It is related to the pionic coupling gB*Bn or 
go'DTT, up to power corrections. The latter is a physically observable quantity that can 
be measured through the strong decay of D* — )■ Dn, and has been discussed at length for 
many years. It also controls the magnitude of loop effects in chiral perturbation theory 
for heavy-light mesons [35] . 

These densities computed on the lattice are presented together with the quark model 
results in Figs. EB and in Tables ITM 

3.2 Difficulties in the comparison between quark models and lat- 
tice QCD 

3.2.1 General remarks 

One is faced with the problem of comparing field theory with a naive quantum mechanical 
model. These are completely different approaches, and it could not be expected that they 
lead to so similar results. Agreement is, however, possible probably because the quark 
models take efficiently into account nonperturbative physical features, like the confine- 
ment, or the spontaneous generation of quark constituent masses. Those features make 
QCD quite different from what one imagines a field theory to be. Indeed field theory is 
often envisaged through the familiar, but inadequate, perturbative picture. However, those 
nonperturbative features can be captured by quark models only in an intuitive way, or in 
average, in an approximate picture. As for perturbative effects, they too can be included 
in quark models (e.g. through radiative corrections, see [36]). It must not be forgotten 
that part of all those effects may happen to be taken into account by the adjustment of 
parameters of the quark models on the experimental data. On the other hand lattice QCD 
(LQCD) is only an approximation to QCD, with its own errors, and with parameters ad- 
justed to reproduce some experimental data. Yet, the essential difference with respect to 

-^^See a discussion in ref. [T7] . 




(43) 



18 



quark models is that LQCD is a systematic approach based on field theory, allowing one to 
approach QCD as much as one would like by increasing the physical volume of the lattices, 
the bare coupling constant, and including all the physical quark fiavors in the background 
gauge field configurations. 

Remark on the scalar density: Physically, the scalar density multiplied by a light 
quark mass perturbation results in a small perturbation of the energy levels of the heavy- 
light state, i.e. 

dE = 5mg X gs . (44) 

In contrast to vector and axial density the scalar one corresponds to a divergent operator 
so that its renormalization constant, Zs, is scheme and scale dependent. Since the mass 
renormalization constant cancels against the one related to the scalar charge, there is no 
need to specify the scheme or scale in the above equation. However, for our purpose where 
we consider the scalar density itself the problem of renormalization should be addressed. 
The renormalized scalar density computed on the lattice has no real intrinsic physical 
meaning, which makes problematic a comparison between the quark models and lattice 
results in this case. For that reason it is more judicious to divide the scalar density by the 
scalar charge which is renormalized by the same constant Zs, and therefore the scale and 
scheme dependence cancels out in the ratio. In other words, we expect 

PQujr) _ Piatt (r) 

QM ~ „latt • V^"^/ 
93 3s 



Disconnected diagrams: In the field theory calculation, with Np ^ 0, there exist 
contributions coming from the insertion of the density operator on an isolated quark loop. 
These are called disconnected diagrams, meaning that the loop is not connected to the 
valence quark lines except by the gluon field. While they vanish automatically for the 
isovector operators, they may be important in certain cases with isoscalar operators. Even 
if they can be sizable, such contributions have been ignored in our lattice calculation of 
densities [7]. Concerning a comparison with quark models, neglecting the disconnected 
diagrams is just what we have to do, since of course the quark models do not include these 
effects in the calculation of currents either. We should simply avoid considering currents 
implying such contributions, which we are of course free to do. The exception is the dis- 
cussion of the meson mass differences due to differences of quark masses. 

Choice of the light quark mass: While the heavy quark mass is taken to be infinitely 
large, a more delicate issue is to decide which light quark mass to consider for comparison 
between quark model results and those obtained in LQCD. There are in fact two distinct 
questions: 

(i) How to fix the mass of the light quark? We should fix it both on the lattice and 
in the quark models to correspond to the same physical spectrum. In particular. 
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using Kg = 0.1374 corresponds to the strange heavy-light meson states. We find that 
the same states are described with cq = —0.1 GeV in eq. ffTOl) [Dirac model\, and 
m, = 0.42 GeV [BT model\. 

(ii) For which hght masses the comparison between lattice and quark model results could 
be made optimal? It turns out that both lattice QCD and quark models are more 
predictive when away from the chiral limit. A comparison is better for the strange 
light quark mass, for which the information on physical spectrum is available too. 
Therefore we choose to work with strange heavy-light meson states. 

Concerning the point b), the cost of numerical simulation of QCD on the lattice becomes 
prohibitively large when going closer to the chiral limit, and reaching the physical m, d 
quarks necessitates chiral extrapolations. The strange quark mass is, instead, directly ac- 
cessible from the current lattice studies. On the other hand, the predictions of our quark 
models can be directly formulated at will for any light quark mass. At present the con- 
stituent quark models do not include the effects of chiral loops which might be significant 
at very small quark masses, such as the case of electric radii of mesons and baryons whose 
expected logarithmic rise at small pion masses is not satisfied by the models [T3] . 

Respective ranges of validity of quark models and lattice QCD: Some com- 
ments are in order regarding the respective domains of validity of quark models and lattice 
QCD as a function of r. 

1) At short-distances, one would expect hard gluons and quantum loop corrections to 
play a dominant role. In particular, the interquark potential should be dominated by 
hard gluon exchange. This effect is only mimicked by the Coulomb potential in the 
quark models, which is an effective potential, with a range much beyond the proper 
UV region, as we already discussed in sec. I2.1.1[ In lattice QCD, on the other hand, 
one should be particularly careful when probing the small distance region where the 
discretization errors might become large. 

2) At long-distances, the confining potential is expected to be modified by pair creation: 
the two-body potential should become fiat, because the energy is now entirely con- 
verted into pairs. This fiattening is indeed demonstrated in lattice QCD for Np > 0, 
around r = 1.25 fm, which means r = 6.25 GeV~^, for the QQ potential, by study- 
ing the mixing between QQ (bb) and di-meson (BB) states [T^. One would then 
expect the quark model density distributions to have too steep a decrease at large 
r, because the confining potential is too strong there. In practice we observe that a 
good agreement between quark models and lattice data extends farther than naively 
expected (up to 4 GeV~^ and even beyond). 

Further difficulty might arise from the fact that we use a lattice with a finite size. More- 
over, due to the periodic boundary conditions, only half of its distances can be explored, 
i.e. r < 6 GeV~^. Fortunately all distributions discussed in this paper are fast decreasing 
functions of r and die out before reaching the boundaries, as shown in ref. [7]. The size of 
the available distances explored in that paper were shown to be enough for this study. In 
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other words, 6 GeV~^ in physical units appears to be already far beyond the mean radius 
of the ground state ^ a/5 GeV~^, and beyond the observed range of agreement of the quark 
model and lattice distributions (around 4 GeV~^). 

3.2.2 Errors in quark models 

In the quark model it is difficult to make a priori statements about errors since it is not a well 
defined approximation scheme to QCD. The only real way to gauge the errors is a posteriori, 
through comparison with experimental data, or with more rigorous theoretical approaches. 
From comparison with experiment, we know that we can obtain a good description of 
spectra and of a certain set of matrix elements, but we know that it is not a quantitative 
statement, that the agreement is very dependent on the observable considered, and on 
the particular model. No model of course can claim an accurate prediction of all possible 
observables. 

Let us just underline two important facts : 

a) Not every quark model shows an extensive agreement with experiment, and in gen- 
eral very few models have been tested on a very large set of data. This is an argument in 
favor of using the GI model which shows good agreement with a very large set of hadronic 
masses and decay widths. Concerning the decay widths, however, it should be stressed 
that the basic eigenvalue equation of GI has to be supplemented by additional modeling 
like the elementary quantum emission model, or the quark pair creation model. For the 
purpose of calculating spatial distributions, we insert it into our Bakamjian-Thomas ap- 
proach to current matrix elements, a combination which has been shown to be successful in 
many places. The success of this approach in predicting important features of Isgur-Wise 
functions is an important argument: it tests both the potential model and the treatment 
of hadron motion. 

More generally speaking, one can say that there is growing evidence, after many years, 
that only relativistic models can provide a consistent picture of processes involving light 
quarks. This applies especially to the internal velocities. 

b) As for the spatial distributions of the static heavy-light mesons considered in the 
present paper, we think a priori that a static limit is a favorable place for quark models; 
and also favorable is that we consider only elastic matrix elements, without considering 
transitions involving radial excitations that would be more sensitive to errors. In particular, 
this is the best place for quark models based on the Dirac equation, while, for instance, 
they would not be suited for annihilation processes. The potential that we use in the Dirac 
equation has not been tested on such a large set of data as the GI one, simply because it is 
specific to heavy-light quark systems. But it appears in practice to be close to the GI one. 

A third favorable condition is to choose as light quark the s rather than the m, d ones. 
This allows, as we said, to avoid as much as possible the phenomenon of spontaneous 
symmetry breaking, beyond the scope of most constituent quark models. 

This favorable context being taken into account, it must be said that the degree of 
agreement which we shall display below for both types of models, is surprising in the 
context of quark models, especially for the (1/2)"'^ state. Quark models are often thought 
of as capable of reproducing the rudimentary features of hadrons, and here we observe 
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a very good agreement with lattice QCD for detailed information on light quark current 
densities as functions of the distance from the static heavy quark source. 

3.2.3 Lattice errors 

We must recall the traditional distinction between statistical and systematic errors. Sta- 
tistical errors are easy to compute. They are often the only ones quoted in the tables and 
figures, but they are most often not the only ones, since the systematic errors are difficult 
to evaluate, and difficult to reduce. 

In the present case, the statistical errors are very small, as it can be seen in Figs. [3|^ 
As for the sources of systematic uncertainties, they can be divided in discretization errors 
and finite volume effects. As we already mentioned, the results obtained in ref. [7] are 
obtained using the lattice box of the size L = 6 GeV~^, which appears to be larger than 
the mean radius of the heavy-light meson states discussed in this paper. For that reason, 
and since we are not working with very light quarks, we believe that systematic errors due 
to finite volume are negligible. 

As for discretization errors, in ref. [15] it was demonstrated on the case of hadron masses 
and decay constants that working with the strange quark mass and with (9(a)-improved 
Wilson action, combined with the Iwasaki gauge action, a bulk of discretization errors 
is minimized. It appears that the results for physical quantities computed in ref. [15] at 
P = 2.1 are very close to their continuum limit, and the difference is always below 10 %. We 
can therefore assume that a similar situation holds true in the present study. Apart from 
"intrinsic" discretization errors, there are also those related to the conversion of results 
from the lattice to physical units, such is the case in eq. P2|) . For that it is necessary to fix 
the lattice spacing by comparing one and the same quantity computed on the lattice with 
its physical value. The most convenient way to fix the lattice spacing is via the Sommer 
radius tq [37] that can be accurately computed on the lattice, but the drawback is that 
its physical value is not known. In ref. [I5] it was shown that the conventionally adopted 
tq = 0.5 fm is compatible with what is obtained if the lattice spacing was fixed by using 
the p-meson mass. Note, however, that the determination of lattice spacing is immaterial 
when computing the adimensional couplings S'vvl 5. 

Finally, we should comment on the sea-quark mass dependence. In this paper we make 
a comparison between the distributions obtained within two classes of quark models and 
those obtained from the QCD simulations on the lattice with both the valence and the 
sea quarks close to that of the physical strange quark. A more physical comparison would 
require an extrapolation in the sea quark to the physical light "u/d" quark mass. We tested 
that changing that mass to half of the strange quark mass (corresponding to Hsca = 0.1382) 
does not change our distributions at all. Whether or not that remains the same when 
the sea-quark mass becomes much lighter is not clear. The CP-PACS Collaboration, by 
working with the same Iwasaki gauge action, checked that the masses and decay constants 
for both light-light [15] and for heavy-light mesons [38] do not depend on the variation of 
the sea-quark mass. A similar conclusion for the heavy-light mesons has been reached by 
using the actions allowing one to probe very light sea quark masses, namely, the improved 
staggered quarks [32] and the twisted mass QCD [ID]. We will therefore assume that a 
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dependence on the sea-quark mass is small. 

For the purpose of comparison with quark models, in the following we will assume that 
the systematic errors of the lattice results are small, at the 10% level, and dominated by 
the discretization effects. The conclusions of the present paper should be understood to 
be valid within the lattice systematics. 

3.3 Quantitative comparison of the two quark models and of lat- 
tice QCD 

Results of the two quark models, and those obtained on the lattice are listed in Tables [H 
El El and displayed in Figs. El HI [5l Following the discussion around eq. ( l45l) in the case of 
scalar density we compare the ratios P^s'\r) / g^g~\ 





r = (1/2)- 


r = (1/2)+ 


r [GeV-^] 


Dirac 


BT 


lattice QCD 


Dirac 


BT 


Lattice QCD 


0. 


0.0534 


0.0536 


0.0446(12) 


0.064 


0.0320 


0.0522(44) 


0.5 


0.0445 


0.0453 


0.0345(8) 


0.0462 


0.0284 


0.0427(34) 


1. 


0.0262 


0.0285 


0.0229(6) 


0.0233 


0.0202 


0.0275(24) 


1.5 


0.0149 


0.0148 


0.0135(4) 


0.0122 


0.0124 


0.0149(15) 


2. 


0.00770 


0.00724 


0.00733(24) 


0.00678 


0.00715 


0.0072(11) 


2.5 


0.00394 


0.00349 


0.00375(13) 


0.00383 


0.00397 


0.00328(51) 


3. 


0.00178 


0.00161 


0.00185(7) 


0.00208 


0.00208 


0.00145(28) 


3.5 


0.000724 


0.000724 


0.000895(38) 


0.00105 


0.00104 


0.00065(20) 


4. 


0.000265 


0.000332 


0.000431(18) 


0.000476 


0.000507 


0.00030(11) 


5. 


0.0000258 


0.0000665 


0.000102(5) 


0.0000683 


0.00011 


0.000076(43) 



Table 1: Comparison between the two models and lattice QCD for py and py in physical units [GeV'^] 



The first striking fact is the similarity of the two quark models for the L = 0, j = 1/2 
state, which was rather unexpected considering that they have a very different formal 
structure. On the other hand, the results concerning the L = 1, j = 1/2 agree less in 
the region of very low values of r, but the overall qualitative agreement is still there. 
Moreover, the agreement at intermediate and larger values of r is very good. We also note 
that contrary to nonrelativistic models, the axial and scalar densities of the L = 1 state 
exhibit nodes in both models, as well as in QCD. 
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r = (1/2)- 


r = (1/2)+ 


r [KjC\ \ 


Dirac 


BT 


lattice QCD 


Dirac 


BT 


Lattice QCD 


0. 


0.0530 


0.0536 


0.0436(10) 


0.0640 


0.0320 


0.0445(39) 


0.5 


0.0377 


0.0437 


0.0313(7) 


0.0360 


0.0245 


0.0320(29) 


1. 


0.0188 


0.0250 


0.0180(4) 


0.0105 


0.0110 


0.0135(13) 


1.5 


0.00969 


0.0116 


0.00890(21) 


0.00153 


0.00276 


0.00273(63) 


2. 


0.00481 


0.00512 


0.00406(17) 


-0.000943 


-0.000192 


-0.00046(62) 


2.5 


0.00224 


0.00219 


0.00179(8) 


-0.00110 


-0.000786 


-0.00087(22) 


3. 


0.000965 


0.000877 


0.00080(4) 


-0.000693 


-0.000616 


-0.00128(24) 


3.5 


0.000377 


0.000346 


0.00037(2) 


-0.000333 


-0.000344 


-0.00167(51) 


4. 


0.000133 


0.000140 


0.000182(8) 


-0.000131 


-0.000164 


-0.00142(73) 


5. 


0.0000123 


0.0000199 


0.0000525(28) 


-0.0000119 


-0.0000250 





Table 2: Comparison between the two models and lattice QCD for pA and p\ in physical units [GeV^] 



It is obvious from the tables and plots that the two quark models agree very well with 
results of lattice QCD. More specifically: 

• Both models are able to reproduce quantitatively the three lattice distributions on a 
large range of distances for both L = and L = 1 states, i.e. 1 GeV~^ ^ ''^ ^ 4 GeV"^. 
Let us recall that there is no adjustment of parameters of the quark models on these 
distributions. The parameters were entirely fixed from the fit to the hadronic spectrum, 
as discussed above. We note however that for the L = 1 states, distributions pre- 
dicted by the models considered in this paper somewhat differ between themselves 
and from the lattice data in the region r < 1 GeV~^. For the ground state, instead, 
the agreement is quantitatively good everywhere, i.e. including the small r-region. 

• Both models reproduce remarkable features of the distributions of excited L = 1, j = 
1/2 states. 

a) From Figs. [3] and H] we observe the similarity of the vector charge distributions for 
L = and L = 1 states as obtained by both models and from our lattice data. Even 
though the overall charge is conserved (fixed to 1), the fact that their distributions 
are very similar was not a priori obvious. 

b) On the contrary, the axial and scalar distributions for L = 1 states are very dif- 
ferent from the ones obtained for L = states. In particular, a change of sign is 
observed for the excited states, a feature not present in the nonrelativistic wave func- 
tions and distributions. Moreover, the positions of zeros are quantitatively predicted 
by the quark models: in the Dirac [BT] model, it is r = 1.7 GeV~^ [r = 1.937 GeV~^] 
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r = (1/2)- 


r = (1/2)+ 


r [KjC\ \ 


Dirac 


BT 


lattice QCD 


Dirac 


BT 


Lattice QCD 


0. 


0.1200 


0.0957 


0.0854(25) 


-0.1251 


-0.1067 


-0.0739(97) 


0.5 


0.0780 


0. 0764 


0.0635(18) 


-0.0608 


-0.0753 


-0.0444(69) 


1. 


0.0341 


0. 0414 


0.0363(11) 


-0.00804 


-0.0217 


-0.0080(34) 


1.5 


0.0161 


0.0179 


0.0174(7) 


0.00745 


0.0069 


0.0090(27) 


2. 


0.0073 


0.0072 


0.00742(48) 


0.00941 


0.0128 


0.0109(23) 


2.5 


0.0032 


0.0028 


0.00298(25) 


0.00698 


0.0105 


0.0073(19) 


3. 


0.0013 


0.00091 


0.00117(12) 


0.00408 


0.0065 


0.0035(11) 


3.5 


0.00045 


0.00028 


0.000456(71) 


0.00200 


0.0034 


0.0015(10) 


4. 


0.00015 


0.000080 


0.000177(33) 


0.000314 


0.0017 


0.00054(30) 


5. 


0.000018 


0.000006 


0.0000246(95) 


0.000102 


0.00031 


0.000050(54) 



Table 3: Comparison between the two models and lattice QCD for ps/gs aiid p^/ds ^'^ physical units 

[GcV^] 

for the axial density, to be compared with r = 1.85 GeV~^ as obtained on the lattice. 
The zero is notably lower for the scalar density, namely in the Dirac (BT) model, we 
find it to be r ~ 1.2 GeV~^ (r = 1.319 GeV~^), once again in a very good agreement 
with the lattice QCD result, r ~ 1.1 GeV~^. Furthermore, from Fig. [5] we observe 
that the signs predicted by our models agree with lattice QCD: in the inner region 
where the wave function is large, the axial (scalar) distribution is positive (negative), 
which is the opposite of what a nonrelativistic calculation would give. 

3.3.1 Explanation of signs and zeroes. Inner versus outer regions 

To explain how these sign changes occur we consider the simpler Dirac formalism. The 
three densities are expressed in terms of upper (/) and lower (g) components of the radial 
wave functions, as given in eqs. (jl]) and ([6]), for L = and L = 1 states, respectively (both 
being j = 1/2). In the nonrelativistic limit , |fi'jy2('")P ^ l/i/2('^)p! which leads to 

Pv{r) = ps{r) = pA{r) = \f[/l{r)\'^, 

(46) 

pUr)=pUr) = \fl/2{r)\^ p\{r) = -l\fl/,{r)\\ 

Pvij) and py{r) are of course positive and, when integrated over the three spatial dimen- 
sions, should give unity. In the nonrelativistic model the same appears to hold true for the 
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scalar densities. The axial densities differ instead in sign (positive for the ground states, 
and negative for the excited ones) and the associated charges also differ in their absolute 
values: for the ground states it is again equal to one, whereas for the excited ones it gives 
— 1/3. The sign of distributions for the ground states remains as such after the relativistic 
corrections are included because \fy\{r)\'^ remains larger than |5']~/2('^)P for all values of r. 
For the excited states, instead, the relativistic corrections drastically modify the form of 
axial and scalar distributions. In the inner region, where the wave functions are most im- 
portant, the relativistic corrections reverse the situation and the upper component |/i/2('^)l 
becomes smaller than |5'J/2('")I- Therefore the axial density is positive instead of negative, 
and the scalar density is negative instead of positive. For small r we note that this dom- 
inance of the lower components is in agreement with the expected behavior, /i/2('") "^^ 
fi'i/2 ^ '"'^ [IS]- For growing r, the upper component |/i/2('")P wins rather rapidly 
over |^!/2(r)P, so that |^i/2(r)P - l/Slfl^^ir)]'^, and \fl/2{r)\'^ - \gl/2{r)\'^, go through zero 
toward the nonrelativistic sign. 

The fact that the zero of the distributions p^^ (r) is closer to the origin than in the case 
of (r) can also be explained: the coefficient of the upper component is larger in the 
scalar case, so the return to the nonrelativistic sign happens earlier. 

3.3.2 A 75-symmetry? The true and the false 

In the Dirac model, in the inner region we observe an exchange of magnitudes between the 
upper and lower components of the radial wave functions in the excited state, with respect 
to the situation in the ground state. Roughly speaking, 

1/1/2(01' ^ 1^2(^)1'' I^V^WI' ^ \fiM^)\"- (47) 

In fact, there is an approximate 75-symmetry, which can be (with suitable normalization) 
written as: 

^(1/2)+ - 75^(1/2)- • (48) 

This symmetry is sufficiently good to explain that the vector and axial distributions of 
(l/2)~ and (1/2)"'" states are similar in this inner region. Lower components give a large 
nonzero value at the origin for all three densities, with the relation dictated by 75 symmetry 
and the vanishing of /i/2(^) at the origin, 

P^(0)=p1(0) = -p+(0). (49) 

These relations, however, are only valid inside what we call the inner region. The situation 
changes at larger r (outer region) where the 75-symmetry does not hold. This may explain 
the failure of the prediction made in ref. [H] concerning the axial couplings, which is solely 
based on chiral symmetry (see below). 

Moments of pA and ps, instead, receive sizable contribution from the large r region 
(outer region), where the sign is opposite to the low r region. The consequence is that for 
axial and scalar density moments, the approximate 75-symmetry is completely broken, as 
noticed first in ref. [B] for the integral of the axial charge. 



26 



The fact that all these highly non trivial features, predicted in a simple manner by the 
relativistic quark models, are verified by lattice QCD seems a big success for quark models. 

We stress the peculiarity of L = 1 states at not too large r, where the ratio between the 
upper and lower components of Dirac spinors is reversed with respect to the nonrelativistic 
limit. Relativity is no longer a correction; it leads to a qualitative change. 

Within the Dirac formalism, we can even interpret this difference between an inner and 
an outer region. The 75-symmetry corresponds to a vectorlike Hamiltonian. The Dirac 
kinetic energy and the "short-distance" Coulomb-like potential are vectorlike, while the 
sum of the confining potential and the mass term is scalar, 

S{r) = ar + c + m . (50) 

Note, however, that the "short" distance is actually not so short, since the strength of 
the Coulomb-like term is rather large, k = 0.75. On the scalarlike side, Cq = m + c < 
counterbalances the positive ar. Quantitatively, we find that in the inner region there is 
an important cancellation in the scalar part S{r), so that the vector one is much larger. 
For example, 

at r = 0.5 GeV"^ : S{r) = 0.025 GeV, V{r) = -1.5 GeV. (51) 

At r = 1 GeV~^, the vector part is still much larger, and around r = 2 GeV~^ they 
are comparable in size. This implies that there is a rather extended region, say for r < 
1.5 GeV~^, in which the Hamiltonian is mostly vectorlike and where the approximate 75- 
symmetry is valid. This in fact is the inner region discussed above. This observation 
may suggest that the Dirac spinors themselves obey this symmetry in the same region. 
For sufficiently large r, however, the scalar potential wins and we return to a situation 
in which 75-symmetry is badly broken and both (1/2)^ and (1/2)+ states may retain an 
approximate nonrelativistic structure. 

3.4 Comparisons of integrated quantities 

We now consider the charges and radii of various distributions O G {V, A, S}, i.e., 

g^+^ = j d^r p^$\r) , and (r^)^ = j d^fr^p^+\r) , (52) 

where, again, the superscript -|- is used when the external states are orbitally excited heavy- 
light mesons = (1/2)+ [L = 1). After the discussion made in the previous subsection, 
it is easy to anticipate as good an agreement between the integrated quantities computed 
on the lattice and within the two quark model frameworks considered in this paper. Their 
comparison, however, bears its own interest: (1) the integrated quantities help to sum- 
marize and quantify the comparison by a few illustrative numbers; (2) these quantities 
are physically interesting, and some of them are even related to physical quantities, such 
as the couplings describing the strong decay to a soft pion. Radii (and higher moments) 
of distributions also may have the interest of reflecting the respective extension of the 
distributions. 
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We note that such quantities are hardly sensitive to the wave function close to the 
origin, although that region has its specific interest as well. In the case of the (1/2)"'' state, 
two distributions change sign, resulting in small values for axial and scalar charges, and 
this is the reason why in eq. ( l52l) we do not define the radii by normalizing the moments 
to their respective charges, as we did in our lattice paper [7]. Furthermore, concerning 
the scalar density, as we discussed in a previous section, in QCD it corresponds to a 
logarithmically divergent quantity that should be properly renormalized. A conventionally 
adopted renormalization scale when discussing the physics of light quarks is /i = 2 GeV, 
that in our case coincides with the inverse lattice spacing. Therefore, as in our lattice 
paper [7] the lattice QCD results for g^^^ refer to MS renormalization scheme, and the 
scale /i = 2 GeV. That, of course, is not an issue in quark models. The results are 
summarized in Tab. |H 





r = (1/2)- 


r - (1/2)+ 


Quantity 


Dirac 


BT 


LQCD 


Dirac 


BT 


LQCD 


(r2)S;) [GeV-2] 


4.82 


5.17 


6.10(8) 


5.96 


5.42 


4.72(42) 




0.63 


0.71 


0.58(2) 


^ 


0.13 


0.08(3) 


(r2)i+) [GeV-^] 


2.72 


2.87 


2.99(9) 


-1.25 


-0.015 


-0.63(37) 




0.44 


0.56 


0.61(3) 


0.51 


0.30 


0.45(14) 


(r^)?^ [GeV-^] 


1.67 


1.72 


2.33(15) 


4.85 


2.73 


3.8(1.4) 



Table 4: Results for moments of the densities discussed in this paper, charges and mean square radii, as 
obtained by using two quark model frameworks, as well as lattice QCD. Note that the vector charges are 
gv = gy = by definition, while the scalar charges and mean square radii inferred from lattice QCD are 
renormalized at ji ~2 GeV in the MS renormalization scheme. 

Radii of the vector distributions : For the charge distribution, the integral should be 1 
on both sides, so there is no check of the quark model by lattice QCD. As for the positive 
definite mean square radii, we see that the results are rather similar for the two states. 
This is in contrast to what one would expect on the basis of a nonrelativistic model, in 
which the centrifugal potential [oc /(/ + l)/r^] suggests that L = 1 states should be more 
extended than L = ones. Once more, this emphasizes the highly relativistic character of 
the dynamics of the light quark. This is particularly the case with the lattice QCD data 
which indicate the L = 1 states being even less extended than the L = ones. 

Integral of the axial distributions: In this case, on the contrary, the integral is physi- 
cally meaningful. For the ground state qa, commonly denoted as 'g, is an analog to of 
the nucleon axial current. ^ has been calculated in a number of ways. In quark models 
its value is known to be lower than its nonrelativistic limit, ^ = 1, due to relativistic cor- 
rections in Dirac spinors, as first stressed by P.N. Bogolioubov [12] • Indeed from our two 
frameworks, and with potentials specified above, we obtain: ^ = 0.63 (Dirac) and ^ = 0.71 
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(BT). These numbers are in very good agreement with what has been obtained after 
integrating the radial distributions computed on the lattice (see also refs. [13]), and are 
much larger than the results inferred from light cone QCD sum rules |44], but compatible 
with the experimentally extracted ^ = 0.59(7) [15] . 

As for the (1/2)"'" state, the coupling is often referred to as g, and it appears to be 
much smaller and of the opposite sign than its nonrelativistic value, g = —1/3, as can be 
seen from Table HI Such a small value is understandable because of the possible cancellation 
between the region in which p\{r) is positive (small r) and negative (large r). Let us stress 
that this cancellation means that the approximate "chiral" symmetry between the (1/2)"*" 
and the (1/2)" states, which is seen for the local distributions in the inner region and 
for the vector radii, is completely broken in the axial case, i.e. ^ <^'g. This cancellation 
is one more striking success of the quark models, already stressed in ref. [6], and it is in 
contradiction with the main assumption made in ref. [H] , that 7j = 7]. 

Radii of the axial distributions : The remarkable fact is that, within the convention 
adopted above flS2]) . the axial radius of L = 1 state is neatly negative both in quark models 
and in QCD. This is different from the sign of Tj because an extra factor inside the 
integral exacerbates the importance of the region of large r, in which p\{r) is negative. 

Integrals and radii of the scalar distributions: As it has been explained above, these are 
not physically meaningful quantities by themselves. Only the ratios, like {r^)s/gs or g^ / gs 
in which the renormalization arbitrariness is eliminated, computed on the lattice and by 
using various quark models, can be compared. In fact, on the lattice we obtained (r^) si 9s = 
3.8(2) GeV~^, that is in very good agreement with what we obtain with the Dirac model, 
{r'^)s/gs = 3.8 GeV~^, and with the BT model, {r'^)s/ Qs = 3.1 GeV~^. Similar agreement 
is observed for the L = 1 state, namely the lattice value {r'^)s/9s ~ 8-4(3.1) GeV~^, agrees 
with the quark model results, (r2)+/^+ = 9.5 GeV"^ [Dirac], 9.1 GeV'^ [BT]. 

We note however that other ratios are in less good agreement. For example on the 
lattice we obtain, g^/gs = 0.74(23), while the quark models give gs/9s = 1-16 [Dirac], 
0.54 [BT]. Similarly {r^)t/{r^)s = 1-6(6), 2.9, 1.6, for the lattice QCD, Dirac, and BT 
models, respectively. 

We note that in contrast with the axial counterpart, the integral gs for the (1/2)"^ state 
is not especially small, in spite of the fact that there is also a zero in Pg{r). It shows 
that the complete cancellation found in the axial case is something nontrivial, resting 
on a quantitative agreement of the structure of the wave function with what is observed 
on the lattice, as it can be seen from eq. (jH]). With the larger coefficient of the upper 

I 2 

components (oc |/i/2('")| ) scalar density the zero appears earlier than in the axial 

case and therefore the small r region [in which Pg{r) is negative] no longer cancels p^ir) 
in the large r region, where it is positive. From the above discussion we see that the sign 
and magnitudes of the integrated quantities for L = 1 states constitute a sensitive test of 
the internal structure of quark models because they probe the lower component of Dirac 
spinors. Moreover, it can be easily understood that (r^)^ is positive, since the positive 

^^For a detailed discussion related to the quark model calculation of g, please see ref. [17]. 

^^While making this comment we have in mind the Dirac model. The same argument also holds also for 
the BT model. Indeed, with a little extra work it can be shown that the BT expressions can be cast into 
the forms similar to the ones written in eqs. (HHl]). 
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contribution of large r is even more favored. 



4 Conclusion 

In this paper we made a quantitative comparison of various distributions of the hght 
valence quark as a function of its distance from the static heavy quark within the heavy- 
light mesons, as previously obtained by using lattice QCD [7], and the corresponding 
quantities in two types of quark models. For this purpose we use a model based on the 
Dirac equation, and another one formulated in the relativistic framework of Bakamjian- 
Thomas. In both models the parameters in the potentials have been fixed to describe the 
hadronic mass spectrum only. As a result we find that the detailed comparison between 
the distributions for L = state [(1/2)^] shows excellent agreement of both models with 
the lattice counterparts. 

Moreover, the agreement also holds for the orbitally exited states, L = 1 [(1/2)"'"], which 
is highly non-trivial because the shape of the distributions are completely different, and 
because the expressions derived within quark models are very sensitive to the relativistic 
corrections and the structure of the potential. 

We focused on three bilinear light quark operators, vector, axial, and scalar, and for 
all of them and for both heavy-light mesons [(1/2)^] we find the agreement to extend over 
a large range of distances, t < 4 GeV~^, in which the distributions of these densities are 
strongly varying (they change values by 2 orders of magnitude). 

What conclusions can be drawn from this agreement? 

- We do not share the widespread opinion that one can get anything at will out of quark 
models because of an alleged large freedom. We repeat that our quark models were 
devised a long time ago, and that no adjustment has been performed on the lattice 
densities which are the object of the present study. In fact there is not much freedom 
in fixing the models if one wants to describe a sufficiently large set of experimental 
data, for instance hadron masses. Moreover, in describing the light quark dynamics 
the inclusion of a relativistic kinetic energy is necessary. It has a clear impact on the 
structure of wave functions. Indeed, the matrix elements of current densities would 
not agree at all with lattice QCD if we used the nonrelativistic treatment, even if 
the spectra could be reasonably well fitted. In particular, one would not get zeroes 
for the distributions of the axial and scalar densities in the case of L = 1 states. 
In addition, the specific Lorentz structure of the potential, scalar at long-distances 
(confining part) and vector at short-distances, seems the simplest way to get the 
small spin-orbit splittings that are required for the static spectrum of heavy-light 
mesons. This structure is particularly important for understanding the behavior of 
distributions considered in this paper in the so-called inner region (r < 1.5 GeV""^). 

- Having adopted the above arguments in constructing the quark models, it may not 
appear too surprising that the resulting radial distributions p^Q\r) {O = V, A, S) are 
similar. However it remains surprising that they agree so well and over such a large 
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range of distances r. Still more surprising is their agreement with the field theoretical 
approach to strong interactions captured by numerical simulations of QCD on the 
lattice. Not only do the distributions have similar shape but they also agree quantita- 
tively on the large range of values of r. In particular, their integrated characteristics, 
(7^^ and (r^)^'', are in good agreement with lattice QCD. 

- Finally, a comment is in order concerning the success of the BT model used in this 
paper. The potential model GI, used in our BT approach, is the one formulated in 
ref. [H]. Its success in describing the detailed lattice data suggests that the wave 
functions derived from the GI model are physically adequate. Therefore, it seems 
likely that their use is not the source of the difficulties encountered in certain quark 
model calculations of heavy to light processes (inclusive rates and form factors) [46], 
ST]. 

We believe that the observations made in this paper further emphasize the continuing 
relevance of quark models for the more intuitive physical interpretation of lattice data. In 
addition, after testing the quark models in a detailed description of the structure of the 
lowest states one can more reliably use them for understanding the experimental results 
which cannot yet be treated by means of lattice QCD. Conversely, lattice data are precious 
since they allow one to test models in situations which could not be accessed experimentally. 
The static case of QCD is a particular example where the lattice data can provide powerful 
tools for discrimination of various models. 
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Appendix 



Dirac Model 

In computations with the Dirac model we choose to work with the orthonormal basis of 
pseudo-Coulombic functions, 



/33/2e-^^41(2/3r). (53) 



(n + 1)! 

To enable the reader with sufficient information necessary to reproduce (approximately) 
all numerical results presented in this paper we also give the expansion coefficients in 

n n 

/'(^) = E ' a'.ir) = ^ c'fU^{r) , (54) 
1=1 1=1 

in the case of n = 15, with /3 = 0.75 GeV, namely, 
f = (1/2)- : 

cP e {0.842464, 0.0518048, -0.0876299, -0.0303635, 0.0110048, 

0.0154663, 0.00653305, -0.000644675, -0.00299685, -0.00252853, 
-0.00156505, -0.00125756, -0.00165585, -0.00213105, -0.00180187} , 

c[^^ E {-0.514038,0.0450279,0.103569,0.0373234,-0.000273779, 
-0.00435446, 0.002004, 0.00634121, 0.00680183, 0.00527419, 
0.00342261, 0.00181708, 0.000444552, -0.000628798, -0.00097132} (55) 

/ = (l/2)+ : 

cP e {0.82558, -0.196543, -0.184425, -0.022037, 0.0312719, 

0.0163965, -0.00538404, -0.0133003, -0.0106817, -0.00538832, 
-0.0020597, -0.00168964, -0.00317236, -0.00454244, -0.00383235} 
G {0.376853, 0.307429, 0.061718, -0.0414061, -0.0369034, 

-0.00895755, 0.00599181, 0.00698936, 0.00275505, -0.000804702, 
-0.00233478, -0.00283731, -0.00339297, -0.00398404, -0.00340649} . 

We stress again that the numerical results presented in this paper are however obtained in 
the much larger and precise basis (n > 40). 

BT Model 

In computations with the BT model, for the ground states, we computed the ground state 
wave function [j^ = (1/2)^] by expanding in a truncated harmonic oscillator basis. 



71=15 



m = E C.i-iyW'/'r JJ^^lY' ii, ]exp{-^ ], (56) 



i=0 



(2i + 1)1/33/2 ^ 1 2/32 / 1 2/32 
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With the masses m2 = 0.419 GeV for the strange quark, and in the heavy quark hmit 
(mi = lO'^ GeV), we find the following values for the coefficients: 

C, e {0.961012542,0.199752584,0.158164397,7.29950815 X 10"^ 
5.70315719 X 10-^ 3.4423735 x 10-^ 2.76433732 x 10-^ 
1.94679499 x 10~^ 1.55471498 x 10"^ 1.18442774 x 10"^ 
9.54456627 X 10"^ 7.95029942 x 10"^ 5.91434911 x 10~^ 
5.39294863 x 10-^ 3.7691386 x 10-^ 3.7654086 x 10"^}. 

Concerning the excited states, = (1/2)+, we expand in the basis: 

and find the coefficients to be 

A e {0.9684052,0.1689021,0.1572283,6.3688330 X 10-^ 

5.2452430 x 10-^ 2.9747754 x 10~^ 2.3702351 x 10~^ 
1.6382594 x 10^^ 1.2990386 x 10-^ 9.7295996 x 10-^ 
7.6334681 x 10'^ 6.3195196 x 10-^ 4.7800541 x 10-^ 
4.1553467 x 10~^ 3.0716809 x 10-^ 2.9703260 x 10"^}. (58) 
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Figure 3: Comparison of the results for three densities [vector pv{t), axial pA{r), and scalar ps{r)lgs] 
for the lowest doublet of the static heavy-light mesons, = (1/2) ~, as obtained by using two quark 
models discussed in the text, namely a Dirac model (dashed line) and a model a la Bakamjan-Thomas 
(solid line). They are displayed together with results obtained in lattice QCD reported in ref [J. 
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Figure 5: Plots showing the region in which the axial and scalar densities in 6.g.\^change the sign, = 
(1/2)+. Dot-dashed horizontal line, corresponding to zero densities, is added for the reader's convenience. 
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